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We derive a many-site version of the non-Mar kovian time-convolutionless polaron master equation [S. Jang 
et al., J. Chem Phys. 129, 101104 (2008)] to describe electronic excitation dynamics in multichromophoric 
systems. By treating electronic and vibrational degrees of freedom in a combined frame (polaron frame), 
this theory is capable of interpolating between weak and strong exciton-phonon coupling and is able to 
account for initial non-equilibrium bath states and spatially correlated environments. Besides outlining a 
general expression for the expected value of any electronic system observable in the original frame, we also 
discuss implications of the Markovian and Secular approximations highlighting that they need not hold in the 
untransformed frame despite being strictly satisfied in the polaron frame. The key features of the theory are 
illustrated using as an example a four-site subsystem of the Fenna-Mathews-Olson light-harvesting complex. 
For a spectral density including a localised mode, we show that oscillations of site populations may only 
be observed when non-equilibrium bath effects are taken into account. Furthermore, we illustrate how this 
formalism allows us to identify the electronic and vibrational components of the oscillatory dynamics. 



I. INTRODUCTION 

Electronic resonance energy transfer is a widespread 
phenomenon in a variety of systems ranging from 
biomolecular components of the photosynthetic 
machinerjff'H', DNAP and fluorescence-based sen soriP 
to conjugate polymeria, crystal impuritie d 9 * 10 ! and 
quantum dot arrayd^MH Traditionally, electronic 
energy transfer in many of these s ystem s has been 
described with Forster-Dexter theory-SCSI which gives 
account of a Pauli-type dynamics of the probabilities of 
individual chromophores being excited or de-excited, as 
if an excitation wer e "ho pping" around. This theory has 
been widely appliecP^CU ana j s particularly successful in 
describing energy transfer when the electronic coupling 
between chromophores is very weak in comparison to 
their interaction with fast-relaxing vibrational degrees 
of freedom. Outside this limit, the electronic interaction 
becomes significant in comparison to exciton-phonon 
coupling and quantum effects can take place, producing 
for example, coll ective (excitonic) behaviour-^*!, coherent 
exciton dynamic d 20 * 21 ! or interference of energy transfer 
pathways^. Therefore, it is of considerable interest 
to have a theoretical framework encompassing in a 
unified manner the various types of electronic excitation 
dynamics that can be observed in multichormophoric 
assemblies. 

The interest in such a unified framework, and par- 
ticularly in coherence effects, has been further moti- 
vated by recent experimental works probing incisively 
ultrafast excitation dynamics and witnessing coherent 
evolutions of excitonic superpositions in light-harvesting 
antennae^HH A common approach to study the ef- 
fects of coherence in exciton relaxation has been the 
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use of Redfield type master equation d 27 * 28 !, where the 
exciton-phonon coupling is the smallest energy scale in 
the open quantum system and therefore can be treated 
as a perturbation to derive a master equation for the 
electronic degrees of freedom. This approach has been 
used to provide insi ghts i nto the effects of coherence in 
spectroscopic signalsPsEII anc j more recently to under- 
stand the roles of coherence and relaxation in the effi- 
ciency of excitation transfer^H^U. However, many multi- 
chromophoric systems operate in an intermediate regime 
where exciton, vibronic relaxation and exciton-phonon 
coupling energy scales are comparable and hence tradi- 
tional perturbative treatments become inaccurate. This 
has lead to recent investigation of excitatio n dynam- 
ics using non-perturbative approache d 39 * 40 * 42 ^^ and so- 
phisticated stochastic treatments^! of the system-plus- 
bath dynamics. Though accurate for small aggregates, 
non-perturbative calculations become very inefficient as 
the system size increases, or in the case of multiple- 
excitations even if the number of chropmophores is small. 
It is therefore of much relevance to develop modified per- 
turbative methodologies^^ that can provide an appro- 
priate qualitative and quantitative account of dynamics 
in the intermediate regime whilst being computationally 
tractable. 

In this work we generalize to multichromophoric ag- 
gregates the polaron-modified mas ter eq uation formal- 
ism pioneered by Abram and Silbe^ESl anc j studied re- 
cently by Jang et alP^ and Nazir^Sl. Although pertur- 
bative, this approach interpolates between the two lim- 
its of weak and strong exciton-phonon coupling, allow- 
ing for a consistent exploration of the regime where the 
energy scales of electronic coupling and exciton-bath in- 
teraction are comparable. In this formalism, the elec- 
tronic system-plus-phonon bath Hamiltonian is trans- 
formed into a new frame (polaron frame) where electronic 
couplings are renormalized and fluctuate due to the in- 
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teraction with the vibrational modes, which are in turn 
fully displaced due to the interaction with the electronic 
excitation. In this way not only is the effect of the bath 
on the electronic system considered, but also a reciprocal 
effect on the phonon bath is accounted for. Under cer- 
tain conditions, the energy scale of the electronic coupling 
fluctuations induced by the displaced vibrations is small 
in comparison to all other energy scales in the system 
and, therefore, such fluctuations can be treated as a per- 
turbation. Standard projection operator techniques can 
then be used to derive a second-order master equation 
that captures non-Markovian and non-equilibrium bath 
effects in the intermediate regime. In addition to the full 
derivation of the homogenous and inhomogeneous terms 
of the time-convolutionless polaron master equation for 
single-exciation dynamics in a multichromophoric aggre- 
gate, this paper contributes the following results: 

(i) We present a general framework for evaluating all 
electronic system observables in the original un- 
tran sform ed lab frame. This extends previous 
woj-lPHlIO which have only been able to consider the 
dynamics of site populations. Full knowledge of the 
reduced electronic density operator in the original 
frame will allow comparison of the non-Markovian 
theory outlined here with experimental data. 

(ii) We illustrate the versatility of the theory to cap- 
ture the dynamical effects induced by structured 
harmonic environments. Particularly, we use as 
an example a four-site subsystem of the Fenna- 
Matthews-Olsen (FMO) complex and show that 
non-equilibrium bath effects, captured by the in- 
homogeneous term, are crucial for the observation 
of long-lived site population oscillations in the pres- 
ence of localised vibrational modes. In general, the 
theory can capture the interplay between electronic 
interactions and localised vibrations in excitation 
dynamics, allowing us to elucidate the possible elec- 
tronic and vibronic origin of the oscillations. 

(iii) We apply two common approximations to the sys- 
tem dynamics in the polaron frame: the Born- 
Markov approximation and the secular approxima- 
tion. Importantly, we are able to show that approx- 
imations made within the polaron frame do not nec- 
essarily hold within the original untransformed lab 
frame; hence, Markovian dynamics in the polaron 
frame may capture some non-Markovian effects in 
the lab frame. Finally, we show that the theory re- 
duces to Forster and Redfield dynamics in the lim- 
its of very weak electronic coupling and very weak 
exciton-phonon coupling respectively. 

The paper is organised as follows. In section |TT] wc 
begin by introducing the multichromophore excitation 
Hamiltonian and then derive a non-Markovian master 
equation describing excitation dynamics within the po- 
laron frame. This section is concluded with a general 



framework for calculating any expected value of the elec- 
tronic system in the untransformed lab frame. In section 
[TTT1 we apply the theory to study the dynamics in a four- 
site subsystem of the FMO complex that interacts with 
both continuous and localised vibrational modes. Sec- 
tion |IV| discusses implications of the Born-Markov and 
secular approximations in the polaron frame, as well as 
the Forster and Redfield limits of the theory. Finally, in 
section [V] we present some concluding remarks. 

II. MANY-SITE POLARON MASTER EQUATION 
A. Polaron-transformed exciton-phonon Hamiltonian 

In this section, we derive a general master equation de- 
scribing single-excitation dynamics of m coupled (chro- 
mophoric) sites interacting with a common bath of har- 
monic oscillators representing the environment (e.g pro- 
tein and solvent). In order to go beyond the weak system- 
environment coupling limit, we perform a polaron trans- 
formation of the exciton-bath Hamiltonian prior to a per- 
turbative expansion with respect to a re-defined system- 
environment interaction in the transformed frame^HSU. 
The Born-Markov approximation is avoided and hence 
both non-Markovian an d non -equilibrium environmental 
effects are accounted foJ 48 * 49 l 

The Hamiltonian describing the combined electronic 
excitation and harmonic environment in an m-site system 
is {h = 1): 

H = Y e rnPm<7m + Y V ™( a m a n + u t °m) 
m (m,n) 

+ y w k^&k + y a ™ a ™ XX 9k . m& k + 

k m k 

(1) 

Here er+ = |m)(0| corresponds to creation of an excita- 
tion on site m with energy e m , V mn denotes the electronic 
coupling between sites m and n, and the notation (m, n) 
signifies a summation over all m and all n > m. The op- 
erator b k (6k) corresponds to the creation (annihilation) 
operator of the k'th mode of the phonon bath, with fre- 
quency Wk- Finally, <?k,m represents the site-dependent 
coupling of site m to the bath mode k with associated 
spectral density J m (oj) = J2k l5k,m| 2 <5(w - w k ). 

We begin our analysis by moving into the polaron 
frame defined by the transformation H = e s He~ s , where 

S = Em ^m^m Ek(«k,m4 ~ «k,™ 6 k), with a k ,m = 

5k.m/ w k- Within this transformed frame, the Hamilto- 
nian for the single-excitation subspace becomes 

H = Y ^n a m + Y ^kfe^k 
m k 

), (2) 

(m,n) 

where the energy of each site is now shifted by its cor- 
responding site-dependent reorganisation energy, A m = 
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Sk uZ ' sucn = £ m — A m . Here we have also 

introduced the new bath operators 

B mn = e^ 5a *— b t- s <^ bk \ (3) 

where <5ak,mn = «k,m — ak,n depends on the difference 
in bath couplings of sites m and n. We now separate H 
into two parts, a non-interacting system and bath Hamil- 
tonian 

H Q = Y l ' mG ™ a m + Y V mnPmn(^m a n + a t a ™) 
m (m,n) 

+ $>k& k 6k, (4) 

k 

and a system-bath interaction described by 

Hi= Y K,»(Bmn^^+5L^%). (5) 

(m,n) 

In doing so, we have defined bath-induced renormali- 
sation factors f3 mn — (B mn ), and shifted bath opera- 
tors B rnn = B rnn — /3 mn . This splitting ensures both 
that (for /3 mn 0) the coherent transfer dynamics 
generated by the bath-renormalized electronic couplings, 
V nm = V nm f3 nm , is fully accounted for within the system 
Hamiltonian, and that (Hj) = 0. For a harmonic oscilla- 
tor bath in thermal equilibrium, the renormalisation fac- 
tors evaluate to p mn = e"5 coth(^ k /2)|5a k , m „| 2 _ Thc 
above implies that in the transformed frame a localised 
electronic excitation fully displaces each mode of the har- 
monic environment, while the environment both renor- 
malizes (see Eq. Q) and causes fluctuations of the elec- 
tronic couplings (see Eq.([5])). Notice that in the limit 
that all sites couple identically to the common bath we 
have that «k,m = CKk and 5a^ mn — for all m and n. 
Hence, all renormalisation factors j3 mn evaluate to unity, 
while the bath operators B mn evaluate to the null oper- 
ator. Therefore, for fully correlated fluctuations we find, 
unsurprisingly, that the dynamics of a single-excitation 
is fully decoupled from the bath. 

To simplify further analysis, we now move into a basis 
in which Hq is diagonal, denoted as the renormalised ex- 
citon basis and labeled with greek letters throughout the 
paper: Ho\a) — e a \a). We may then express the original 
system operators in this new basis: <r+ = ^ Q u ma a+ 
where a\ = |a)(0| and u ma = (a|r7i). The polaron- 
transformed Hamiltonian within this new basis becomes 

Hq = Y ^ a t a « + Y ^k&k&k, (6) 

a k 

Hi = Y E ( V ™ B 

m n 'Urn a ^ n Q & a & @ 

otfi (m,n) 

)■ (7) 

Moving into the interaction picture defined with respect 
to Hq is straightforward, and we obtain the interaction 



Hamiltonian 

a/3 (m,n) 

+V mn Bl n (t)u^ a u nP a+a-e^-^ t y8) 
= ]T (Sap® ® B a0 (t) + Slp(t) ® Bl fj (t)) , 

a/3 

(9) 

where we have defined new system and bath operators 

B a p{t) = Yj VmnUmaUnpBmnit), (10) 
(m.n) 

with e af} = e a ~ep and B mn (t) = e tflot B mn e- iilot . We 
therefore have a simple form for the polaron-frame in- 
teraction Hamiltonian Hj(t), which will be treated as a 
perturbation in the master equation derivation that fol- 
lows. 

B. Time-Local Master Equation 

The time-convolutionless and the Nakajima-Zwanzig 
master equations are, respectively, time-local and time- 
non-local perturbation expansions to describe non- 
Markovian dynamical evolutions^. However, there are 
several reasons one can chose the time-local expansion 
over its non-local counterpart. As is illlustrated by Vac- 
chini and Breuer^, these two approaches have different 
ranges of validity: while the time-convolutionless expan- 
sion breaks down at finite time in the strong coupling 
limit, the Nakajima-Zwanzig approach does not necesar- 
ily preserve positivity if restricted to second order. This 
suggests that if one is able to identify the parameter space 
of weak system-bath coupling (as we do in our case), 
a second-order time-convolutionless perturbation scheme 
may be appropiate appropriate. Furthermore, it is well 
known that time-local approaches are often simpler. This 
is particularly important in our case given that the po- 
laron treatment for many-sites is already quite involved. 
Most importantly, however, current research is helping 
to clarify that time-nonlocal and time-local approaches to 
non-Markovian dynamics are complementary rather than 
opposed. For instance, Chruscihski and Kossakowski 59 
show the retarded time integration in a time-nonlocal for- 
malism can be mapped onto a local- time dynamics with 
a generator that has a strong dependence on the starting 
point to • Because of the above, in this work we consider 
a time-convolutionless projector operator expansion to 
derive a second-order, time-local master equation in the 
polaron representation. 
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1. Projection Operator Formalism 



2. Homogeneous Term 



Having transformed our Hamiltonian into the polaron 
frame, we now wish to derive a time-local master equa- 
tion governing the reduced dynamics of our multichro- 
mophoric excitation under the influence of the harmonic 
environment. In order to do so, we follow the time-local 
projection operator formalism (as given, for example, in 
Breuer and Petruccione^) . In brief, we define a projec- 
tion super-operator V as 



X -> Vx = tr B {x} ® Prcf 



(11) 



which projects onto the relevant part of the combined 
system-environment density matrix \i such that Vx gives 
the complete information required to reconstruct the re- 
duced density matrix of the open system. Here, p Te { de- 
notes a fixed (arbitrary) reference state of the environ- 
ment, commonly chosen to be the thermal equilibrium 
state. The complementary super-operator Q is also de- 
fined, through Qx — X ~ T^Xi which projects onto the 
irrelevant part of the density matrix. 

Suppose the dynamics of the combined system plus 
bath is governed by a Hamiltonian of the general form 



H = H + aH I , 



(12) 



where Hq determines the uncoupled time evolution of 
the system and environment, Hj describes the system- 
environment interactions, and a denotes a dimensionless 
expansion parameter. By applying the projection opera- 
tors to the interaction-picture Liouville equation 



d X (t) 
dt 



-ia[H I {t),x{t)]=a£(t) X {t), 



(13) 



we may derive an exact time-convolutionlcss master 
equation for the relevant part of the density matrix of 
the form 



dt 



V X (t) =K(t)x(t)+I(t)x(t ), 



(14) 



provided neither a nor t — to become too large, where to is 
the initial time. To second order in the expansion param- 
eter a, the homogeneous term lZ(t) and inhomogeneous 
term I(t) can be given as 

K{t) = a 2 [ ds VC{t)C{s)V, (15) 
Jo 

l(t) = aV£(t)Q + a 2 [ ds PC(t)C(s)Q, (16) 
J o 

respectively, where we set to = and have imposed the 
condition that VC{t)V = 0. Here, the inhomogeneous 
term I(t) is non-zero only when the initial bath state 
differs from the reference state p Te f . In our case, we shall 
see that this corresponds to a non- equilibrium prepara- 
tion of the initial environmental state within the polaron 
frame. 



Let us now look at the homogeneous term of our mas- 
ter eqaution for the polaron frame interaction Hamilto- 
nian [Eq. ph]. To obtain the reduced dynamics for the 
many-site system, p(t) , we trace over the bath degrees of 
freedom to arrive at: 



K(t)p(t) =tr B {a 2 



ds VC(t)C(s)Vx} 

ds tr B {[C(t)C(s)p(t)Prc t ]} 

ds tY B {[H I (t),[H I (s),p{t)p Icf }]Y 

(17) 

We shall now take our bath reference state p rc { to be the 
thermal equilibrium state within the polaron frame i.e. 
p re f = pb- On substituting in the interaction Hamilto- 
nian given in Eq. d9]), we obtain 

R(t)p(t) = - £ (^Lu(t)e i£ ^[Sap,S^p(t)] 
+^,^t)e^ t [S aP ,Sj af p(t)] 

+^l^ty^ t [su,sij(t)] 

+h.c.Y (18) 



Here, we have defined the time-dependent rates 



(i) 



Jo 

JO 

Jo 

with corresponding two-time bath correlation functions 

c£U('-«) = £ (B mn (t)B pq ( S )), 

{mn) {pq) 

c aLuit-s) = £ y.k;p (BUt)B pq ( S )), 

{mn) {pq) 

<$L(*-") = E H u ™ (B mn m q (s)), 

{mn) (pq) 

C tL(t- s ) = £ H U ™ (B-Ut)Bi q (s)), (20) 

{mn) (pq) 

where U™™ = V mn V pq u ma u n pU w u qv , and (...) denotes 
the average with respect to the reference state of the bath 



(2) 



(t- 


*), 


(*- 


s), 


(*- 


s), 


(t- 
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p B . The derivation of the explicit form of the two-time 
correlation functions is rather involved, so we show the 
exact expressions in Appendix |A"| 



3. Inhomogeneous term 

If the initial bath state within the polaron frame differs 
from the reference bath state p B i-e. it is a non-thermal 
equilibrium state within the polaron frame, then we must 
also account for the inhomogeneous term in the master 
equation (see Eq. (14)). 

Let us consider a general separable initial state in 
the lab frame (i.e. prior to polaron tranformation) : 
X(0) = Pij(ty a t a 7 ®Pb, where p B denotes the ther- 
mal equilibrium bath state in the lab frame. Transform- 
ing into the polaron frame we find the initial state 

X(0) = J^p tl (0)<7+ajl[(3- J 1 D(a k , l )p B D(~a Kj ). 

ij k 

(21) 

Here Pij(0) = PijPij(0) and denotes the ij'th element of 
the initial system density operator in the polaron frame 

is the bath displacement 



l 2 (t)p(0) =a 2 ds tr B {P£{t)£(s)Qx(0)} 
Jo 

= - f dstrsl^), [#/(*), Qx(0)]]}, 



(25) 

Substituting in the interaction Hamiltonian, we obtain: 
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PiM^l^uity^iS^, S^afa 



3 ' 

(26) 



a kJ bl-a 



and D(pik j) = e 
operator of mode k due to interaction with site j. 

The irrelevant part of the total system-bath density 
matrix at time zero is then given by 

Qx(0) = X>i(0)o-+*7 II (&D{a^~p B D (-okj) - p B 

ij k 

= Y t PiM°t a iQiiPB- (22) 

ij 

Notice that we have defined Qijp B as the state account- 
ing for the difference between the displaced bath and the 
bath thermal equilibrium in the polaron frame. We are 
now in a position to evaluate the inhomogeneous term; 
tracing over the bath we obtain 

ds tT B {VC(t)jC(s)Qx(0)}. (23) 



where we have defined the rates: 

fl (*)= / dse^DY-) R (t-s), 







" (2) (t) 
"(3) {f \ 

" (4) (t) 



ds e 



= / ds 



ij,aP,iu 
ij,a 

Mvu 1 n( 3 ) 



(*- 


*)> 


(t- 


s), 


(*- 


s), 


(*- 


s), (27) 



with bath correlation functions 



D 



(i) 

ij,ap,^ 



D 



(2) 



Let us consider each term separately. On substituting in 
the form of the interaction Hamiltonian to the term first 
order in a, we obtain 

M^piO) = -itv B {[Hj(t),Qx(0)]} 



„(*-*) = E Y. U "{B mn (t)B pq {s)) QirPB , 

(H <P<?) 

(t - «) = E T,K^(Bln(t)B pq {s)) Q%] ~ PB1 

{mn) {pq) 
{mn) {pq) 

„(*-«) = E T,^{Bl n {t)Bl q {s)) Qij ~ PB . 



D 



(4) 



{mn) {pq) 



(28) 



+ -71 
3 1 



ap ij 

+ h.c 

Here we have introduced the rate T 



(24) 



S(m,n> VmnUmaU* n p{B rnn \T, ) ) Q tj p 



(*)>« 



- 7, a/3 (^) 

see Appendix |B|, 



where (■■■)o, ij p B denotes a thermal average with respect 
to Qijp B . 

The second order in a term of the inhomogeneous 
super-operator is 



Once again, explicit forms for these correlation functions 
are presented in Appendix [B] 



C. Lab Frame Dynamics 

The master equation derived above gives the dynamics 
of the reduced density matrix for the electronic system 
within the polaron frame. However, we are interested 
in the excitation dynamics in the original untransformed 
lab frame. To calculate the correct transformation from 
polaron to lab frame, consider the Schrodinger picture 
system-bath density operator in the polaron frame x(t) = 



G 



e s x(t)e ~ s . Inverting this expression and using the iden- 
tity V + Q = I, we may write the lab frame combined 
density operator as: 



x (t) = e- s Vx(t)e S + e~ s Qx(ty 



(29) 



The expectation value of a system observable A in the 
lab frame is given by 



(A) = tv s+B {A x (t)} 

= tr s+B {e s Ae- s Vx(t)} + tr s +B {e s Ae~ s Qx{t)} 
= (A) Iel + (A) fed . (30) 

From the definition of the projection operator, the first 
term (A) re i is trivial to evaluate: 

(A) lel = ti s +B{e s Ae' s Vx(t)} 

= ti s+B {e s Ae~ s p(t) ® p B } 

= tt s {Ap(t)h (31) 

where we have defined the transformed observable A — 
tr B {e Ae~ p B }. Since this contribution depends en- 
tirely on the relevant dynamics, we have defined it as 
the relevant contribution to the expected value of A. 

To evaluate the second term, (A}- lITe i, we require knowl- 
edge of the dynamics of the irrelevant part of the den- 
sity matrix. Breuer and Petruccione^ show that the 
irrelevant part at an arbitrary time t can in principle 
be determined from the knowledge of both the rele- 
vant part Vx(t) and the initial condition Qx(0). Im- 
portantly, in our case, both these quantities are known. 
Therefore, we may formally write the irrelevant part as 
Qx(t) = S(t)Vx(t)+T(t)Qx(0)- The resulting irrelevant 
contribution to the system operator expectation value is 



(A) irre] = tr s+B {e s Ae- s Qx(t)} (32) 
= tr s+B {e s Ae- s S(t)Tx(t)} 
+tT S+B {e s Ae- s T(t)Qx(0)}- 

(33) 

Up to second order in the coupling a, the super-operators 
S(t) and T(t) are given by: 

S(t)=af dsC(s)+a 2 f dsf ds 1 QC{s)C{s'), 
Jo Jo Jo 

T(t) = l + a[ dsQ£(s) + a 2 [ dsf ds' QC{s)QC{s') 
Jo Jo Jo 

-a 2 f ds f ds'Q£{s')V£(s). (34) 
Jo Jo 

In general, the exact forms and calculations of the irrel- 
evant contributions to expectation values are extremely 



involved and for simplicity here we restict ourselves to 
operators for which such terms evaluate to zero as it is 
explained below. However, it is worth noting that to 
zeroth order in the coupling parameter a the irrelevant 
contribution becomes {A)- lm \ = tvs+ B {e s Ae^ s Qx(0)} , 
which can be used to evaluate approximate expectation 
values of system operators that do not commute with the 
polaron transformation S. Specifically, the expectation 
value of an observable in the lab frame including just the 
zeroth order term for the irrelevant contribution reads 



(A) = tr s {Ap(t)} + tr s {^p(0)} - tr 5 {ip(0)}. (35) 

In the case of system operators commuting with the 
polaron transformation S, such as the to— th site popula- 
tion operator cr+er~ , the irrelevant contribution vanishes 
and the expected value in the lab frame is entirely deter- 
mined by the relevant contribution. In other words, site 
populations remain unaffected during the transformation 
back to the lab frame. Let us demonstrate this explicitly: 



Wm<r m ) = t,[ s+B{o-+a- m Vx(t)} + tr s + B {o-+a m Qx(t)} 
= ^s{^ m p(t)} + tr s{CT+er m tT B {Qx(t)}} 
= tesWUvmPit)} = {o-+o-~) ieh (36) 

where by definition trs{%} = trgjT ^x} and VQ = 0, and 

we show that the 
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therefore (c^cr~)i rro i = 0. In Sec. 
key aspect of the theory can be illustrated by considering 
site populations. 

One may attempt to compute off-diagonal operators in 
the site basis, i.e. cr+cr~ with m ^ n, with the approx- 
imation proposed in Eq.(35) which we present here for 
completeness. 

(°"mO = Pmntrs{o-£ l o--p(t)} + tr s {a+a~p(0)} 

-f3 mn tr s {<T+o--p(0)}. (37) 

Let us finish this section with a discussion of the va- 
lidity of this polaron treatment. As our approach is 
pertubative we expect the master equation to be valid 
only within certain regimes. Firstly, notice that in the 
absence of electronic couplings the interaction Hamilto- 
nian Hi (Eq.([5])) is zero and the polaron transformation 
exactly diagonalises the combined system-bath Hamilto- 
nian. Therefore, we expect this perturbative treatment 
to be a good approximation in the limit where the mag- 
nitude of electronic couplings are small in comparison to 
the detunings between onsite energi es, irr espective of the 
strength of the coupling to the batr l 52 ! 53 !. Moreover, this 
perturbative treatment is valid if the energy scale associ- 
ated to the fluctuations of the electronic couplings repre- 
sented by Hj is the smallest energy scale in the system. 
Such fluctuations are given by^ 



|2\ 1/2 



-V (l-B 2 ) 1/2 



(38) 
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Furthermore, as mentioned before, the small polaron 
transformation S assumes that bath modes are fully dis- 
placed by the interaction with a localised electronic ex- 
citation. For super-Ohmic spectral densities such full 
displacement is only valid for bath frequencies larger 
than the typical energy scale of the renormalized exci- 
tonic HamiltoniarPSMl Notice also that for an Ohmic 
spectral density, i.e. J(w) oc u>, and independent baths 
for each site, the renormalization factors exhibit a well- 
known infra-red divergence^ that leads to electronic cou- 
plings being renormalized to zero independently of the 
strength of the system-environment interaction. Both of 
these short-comings can in principle be alleviated by con- 
sidering a variational-polaron approach^ or by develop- 
ing perturbative treatments based on al ternat ive trans- 
formations of the system-bath interactio n 142 ! 43 !. 



III. NON-MARKOVIAN DYNAMICS 

To illustrate the scope of this theory, we apply it 
to study the dynamics of a subsystem of the Fenna- 
Matthews-Olsen (FMO) complex. In particular, we con- 
sider the subsystem involving sites 1, 2, 3 and 4 with an 
electronic Hamiltonian taken from Cho et al¥^. In units 
of cm -1 this reads: 



H, 



FMO — 



280 -106 8 -5 

-106 420 28 6 

8 28 -62 

-5 6 -62 175 



(39) 



Each site is coupled to an independent bath with spec- 
tral density given by J{uj) — so Jq(oj) + sh Jh (wp^, which 
has a continuous contribution Jq(ui) and a contribution 
from a localised vibrational mode Jff (w). The continuous 
part of the spectral density is defined as: 



J (w) = 



si + s 2 



7!2w 4 ' 

i=1.2 1 



-(ui/u 



(40) 



The localised vibrational mode is commonly described 
by a delta function. This provides a simple picture for 
describing the coupling to localised modes and allows 
for analytical expressions for correlation functions, but 
does not necessarily represent a realistic situation. One 
would expect that in practice the single frequency mode 
is bro adene d by interactions with the surrounding bulk 
m0( j e ^55i56| Therefore, in this work we shall assume a 
broadened vibrational mode with a Lorentzian line shape: 



2uj 



H 



0J 3 € 



7r (up- — io 2 H ) 2 + e 



2,. ,2 ' 



(41) 



Here, the parameters for the continuous part of the 
spectral density are sq = 0.5, si = 0.8, s 2 = 0.5, 



u)\ = 0.0069 meV and uj 2 = 0.024 meV. Meanwhile, 
the parameters for the localised mode are sh = 0.22, 
ujjj = 180 cm -1 and broadening e = 50 cm -1 , a value 
that has been chosen to be larger than the average elec- 
tronic coupling strength. For all the calculations, room 
temperature is assumed, i.e. kgT = 200 cm -1 

For the FMO spectral density introduced above, the 
bath renormalisation factors may be written as 



p mn = exp( - j coth03w/2)) 

J h (oj 



x exp 



dw- 



coth(/3w/2) (42) 



and the renormalised electronic Hamiltonian then evalu- 
ates to 



H 



FMO 



280 
-0.107 
0.008 
V -0.005 



-0.107 
420 
0.028 
0.006 



0.008 -0.005 

0.028 0.006 

-0.062 

-0.062 175 



(43) 



The theory predics that the harmonic bath renormalises 
the electronic couplings to such a degree that the exci- 
tonic eigenstates are effectively localised on sites. Hence 
the transition frequency between the two highest en- 
ergy eigenstates is set by the difference between the 
energies of sites 1 and 2. It is worth noting here 
that these strong bath-induced renormalisations predom- 
inantly arise from the continuous part of the spectral den- 
sity Jo(u>). For the parameters given above onsite energy 
gaps are all larger than electronic couplings between sites 
i.e. \e m — e„| > \V mn \ though electronic coupling fluc- 
tuations "f m n are comparable with V mn . Nevertheless, 
all 7 m n are still smaller than the characteristic frequency 
of Jq(uj) (around 200 cm -1 ), which makes the present 
polaron treatment appropriate. 



A. Non-Equilibrium Effects 

We now investigate the dynamics of the FMO subsys- 
tem as predicted by the non-Markovian polaron theory. 
We begin by considering an initial state localised on site 
1. Figure [T] shows the dynamics of populations of the 
four sites and compare the effects of the different terms 
of the non-Markovian polaron master equation. In the 
presence of the homogeneous superoperator term alone, 
all four populations evolve monotonically with no dis- 
cernable oscillatory dynamics. On including the inhomo- 
geneous superoperator in the polaron master equation, 
we see remarkably the emergence of oscillatory dynamics 
in the populations of sites 1 and 2. This behaviour is long 
lived, lasting upto 600 fs. Beyond the 600 fs timescale 
we see no variation between the dynamics considering 
solely the homogeneous term and the full polaron master 
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FIG. 1. Population dynamics of the FMO subsystem assum- 
ing an initial state localised on site 1. Presented are the dy- 
namics with just the homogeneous superoperator (dashed), 
and the homogeneous plus inhomogeneous superoperators 
(solid). 



FIG. 2. Population dynamics of the FMO subsystem assum- 
ing an initial state consisting of a symmetric superposition 
of sites 1 and 2. Presented are the dynamics with just the 
homogeneous superoperator (dashed), and the homogeneous 
plus full inhomogeneous superoperators (solid). 



equation. Therefore, one can conclude that the inho- 
mogeneous terms, which describe non-equilibrium bath 
effects, have a profound effect at short times allowing for 
the emergence of oscillatory dynamics in agreement with 
previous works^H 



B. Delocalised Initial states 

We now consider an electronic excitation initially delo- 
calised over a number of sites while being separable with 
the thermal equilibrium bath in the lab frame. Upon 
transformation into the polaron frame, this state maps 
onto an initially correlated system-bath state. Assum- 
ing an electronic excitation symmetrically delocalised be- 
tween sites 1 and 2, figure [2] depicts the population dy- 
namics of each site, comparing the evolutions given by 
the homogeneous term versus the full the polaron mas- 
ter equation. In the presence of the homogeneous term 
alone, the populations evolve incoherently as in the case 
of localised excitation. Upon the inclusion of the inho- 
mogeneous terms we once again see the emergence of 
well-defined, long lasting oscillations. Interestingly, along 
with the coherent oscillations in the dynamics of sites 1 



and 2, we are also able to observe subtle oscillatory be- 
haviour in the population of site 3. These oscillations can 
be seen to decay over the same 600 fs timescale observed 
in the dynamics for a localised initial state. Therefore, it 
would appear that, for the parameters given, delocalized 
excitations do not have a profound effect on the timescale 
over which oscillatory dynamics is observed. 



C. Role of the localised vibrational Mode 

In this section we explore the effect of the localised 
mode on the population dynamics and illustrate how 
this formalism allow us to identify the vibronic or elec- 
tronic origin of the observed oscillatory dynamics. Fig- 
ure [3] presents the full non-Markovian dynamics in the 
absence and presence of the broadened localised mode. 
An important characteristic of this mode is that its en- 
ergy matches the energy transitions between the highest 
excitonic eigenstates. This resonant condition leads to 
a dramatic effect on the site population dynamics as it 
not only enhances oscillations in the probabilities of hav- 
ing sites 1 and 2 excited, but also it increases the rate of 
energy transfer to lower energy sites 3 and 4. We have al- 
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FIG. 3. Population dynamics of the FMO subsystem with 
(solid) and without (dashed) the localised vibrational mode. 
The dynamics presented assumes an initial state localised on 
site 1. 



ready seen that when considering the full spectral density, 
including the broadened localised mode, there are strong 
long-lasting oscillations in the populations of sites 1 and 
2. However, if the localised energy mode is neglected, 
such oscillations are not present. This agrees well with 
recent results of Prior et aL 43 who have predicted, us- 
ing time-adaptive density matrix renormalisation meth- 
ods, similar strong enhancement of coherent oscillations 
upon the inclusion of a localised mode. In Figure [4] we 
present the Fourier transform of the population of site 
1. In the presence of the localised mode we clearly see 
a strong peak at approximately 180 cm -1 , correspond- 
ing to the energy of the localised mode. Therefore, we 
may associate the observed oscillations to a vibronic- 
induced effect. Meanwhile, in the absence of the localised 
mode we observe a very broad peak centred at approxi- 
mately 150 cm -1 , corresponding to the energy difference 
between the two highest renormalized electronic eigen- 
states. The broad nature of the peak is associated to 
very-short lived oscillatory dynamics. 




a) (cm 1 ) 



FIG. 4. Fourier transform spectrum of the population dy- 
namics of site 1 in the absence (dashed) and presence (solid) 
of the localised vibrational mode. 



IV. MARKOVIAN AND SECULAR APPROXIMATIONS 
IN THE POLARON FRAME 

In general terms, for open quantum systems, non- 
Markovian behaviour describes a system's dynamics 
that is largely affected by bath memory effects and 
the associated system-bath correlations. Given the 
variety of approaches that one can follow to describe 
such non-Markovian phenomena, in the recent years 
there have been an increased number of works at- 
tempting to formally derive a measure for the degree 
of non-Markovianity of a quantum dynamics^HUl. I u 
particular, as it was discussed in RefpS, while Markovian 
processes tend to continuously reduce the distinguisha- 
bility between the dynamcal evoluiton of two different 
initial states, non-Markovian dynamics exhibit periods 
in which there is a growth of this distinguishability. 
This fundamental difference between Markovian and 
non-Markovian behaviour is illustrated in this section 
where we show that an interesting aspect of the many- 
site polaron theory is that Makovian evolutions in the 
polaron frame can actually capture some non-Markovian 
phenomena when transforming back to the lab frame. 
Similarly, we will discuss that a secular approximation 
in the polaron frame does not nccesarily imply that 
eigenstate populations and coherences in the lab frame 
are decoupled. 



A. Markovian approximation 

In order to consider the Markovian dynamics in the 
polaron frame we begin by assuming the Born approx- 
imation which implies that the system and bath states 
factorise at all times in the transformed frame i.e. x(t) — 
p(t) ® ps ■ Here the bath ps is assumed to be in ther- 
mal equilibrium in the polaron frame at all times. As a 
result, we see straightaway that this leads to the inho- 
mogeneous terms given in Eq.(23) to evaluate to zero. In 



conjuction with the Born approximation, we next assume 
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that the bath relaxes on a time scale much shorter than 
the characteristic timescale of the system's evolution i.e. 
Born-Markov approximation. Therefore, we may extend 
the upper-limit of the integrations in the rates in equa- 
tion ( 19 1 to infinity. Now, the rates still have an explicit 
reference to the starting time t = 0. This dependence on 
the past can be made explicit by making the substitution 
s — > t — s. The resulting Markovian master equation is 
d ^jr- = TZ-M(t)p(t) where the Markovian super-operator 
7?.A/(i) is defined as 



Hum® = - E ^al^^^is^s^m] 



^l^ Mt K,, sum] 

+h.c) . (44) 



The time-independent, Markovian rates are: 



.(i) 

" a/3, {iv 



(2) 

a/3, {iv 
(3) 

a/3 , /if 
" a/3,/i,z/ 



ds 



ds e-^° C^(s), 
ds e-«- s ci%^(s), 
ds e~^' C^a), 



(45) 



This describes the Markovian dynamics of electronic 
degrees of freedom dressed with a phonon reservoir. No- 
tice, however, that separability of the system-plus-bath 
state in the polaron frame does not imply separability in 
the lab frame, as the transformation back to the lab frame 
converts the dressed electronic system into bare elec- 
tronic degrees of freedom correlated to the vibronic de- 
grees of freedom through the appropriate displacements. 
We therefore expect that the Markovian polaron theory 
should be able, under certain conditions, to capture some 
non-Markovian effects in the lab frame. To explore this 
idea, we investigate signatures of non-Markovianity dur- 
ing the excitation dynamics. 

A number of recent papers have proposed several mea- 
sures to quantify the non-Narkovian character of quan- 
tum dynamics by quantifying the deviation of dynamical 
evolution from a Markovian limitP^H. Q ne particularly 
simple measure proposed by Breuer et alP^ is based on 
the distinguishability of states as measured by the trace 
distance^. The trace distance between two states p\ and 
P2 is defined as: 
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FIG. 5. Dynamics of the trace distance measure and its 
derivative using a Markovian polaron master equation. Pre- 
sented are the measures evaluated in the polaron frame (dash- 
dotted) and the lab frame (solid). 



where \A\ — V AfA. It can be shown that all trace 
preserving completely positive maps are contractions of 
D. Therefore, all Markovian processes lead to a con- 
tinuous reduction of the distinguishability between dif- 
ferent states and the loss of distinguishability may be 
interpreted as the irreversible flow of information from 
the system to the environment)^. In contrast, a non- 
Markovian process may be characterised by periods of 
increasing distinguishability during the evolution where 
the rate of change of D(pi(t), P2(t)) acquires positive val- 
ues. 

To illustrate these ideas in a specific example, we con- 
sider the same four-site subsystem of FMO. However, in 
order to allow us to make a valid Markovian approxi- 
mation, we consider only the continuous component of 



D( Pl (t),p 2 {t)) = -tr| Pl (t)-p 2 (t)|, 



(46) 



the spectral density (see Eq. (40)) and ignore the lo- 
calised mode that inherently leads to non-Markovian dy- 
namics. Furthermore, we assume that the characteris- 
tic bath frequencies Wi^ are increased by a factor of 10, 
while simultaneously the weightings si,2 are reduced by 
the same factor. This ensures that the bath correlation 
functions decay rapidly, while still maintaining the same 
reorganisation energy A = j du>Lu^ 1 J (uj) characterising 
the strength of the system-bath interaction. 
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Taking p\ and pi as the states where excitation is 
initially localised on sites 1 or 2 respectively, figure 
[5|a) depicts the trace distance between pi and p 2 
as a function of time, evaluated in both the polaron 
and lab frames. Although oscillations of the trace 
distance are observed within both frames, we notice 
that in the polaron frame the trace distance is always 
decreasing, while in the lab frame there are periods 
where this measure increases. To illustrate this point 
further, figure [5jb) presents the derivative of the trace 
distance measure. Clearly, in the lab frame there are 
intervals with positive derivative denoting an increase in 
distinguishability and thus suggesting a non-Markovian 
process. However, in the polaron frame the derivative is 
always strictly negative confirming the Markovianity of 
the dynamics within this frame. Therefore, as expected, 
Markovian dynamics within the polaron frame can give 
rise to non-Markovian evolution in the original lab frame. 



B. Secular approximation 

When the exponential terms of the form exp[i(a;+u/)t] 

average to zero 



44 



in the superoperator given in Eq. 
on the timescale relevant to bath relaxation processes, 
the secular approximation can be made. In that case, 
only the terms where u) + <J = are retained and the 
relaxation super-operator becomes 

Kswm = -E { T alp a [s a p,s Pa m] 



a/3 



+^aLp[SU,S a pp(t)} 

+rS iQ)8 [s„p,s^p(t)] 
+^lp a [slp,slj(t)] 

+h.c. 



(47) 



An important consequence of the secular approxima- 
tion is that in the polaron frame, eigenstate popula- 
tions and coherences are decoupled from each other. To 
show this consider the evolution of the expectation value 
Pnv(t) = {lAP(f)\ v ) ■> where and \v) denote renormal- 
ized excitonic eigenstates. The population /5 w (i) within 
the polaron frame evolves according to: 



dt 



dt 



h.c. 



fl.flOt 



r(2) 

ap,,afi 



Ft 2 ) 
jLtat,£tot 



r(3) 



a 

+h.c 



, p(2) , p(3) 



(49) 



Notice that eigenstate populations obey a simple 
Pauli master equation. Hence, no oscillatory dynamics 
is expected for renormalized exciton populations in the 
polaron frame. However, since the polaron transfor- 
mation does not commute with the rotation between 
the site basis and the renormalised eigenstate basis, 
there is a mixing of populations and coherences upon 
transforming back into the lab frame, and hence beating 
of the population of renormalised excitonic states in 
the lab frame should be observed. To illustrate this 
effect, we consider the same modified FMO system as 
described in Sec. |IV A| Figure [6] shows the dynamics of 
the renormalized exciton eigenstate populations in both 
polaron and lab frames, as predicted by a polaron master 
equation with secular approximation for an excitation 
initially localised on site 1. While population transfer 
proceeds in an inchoherent manner in the polaron frame, 
in the original lab frame we observe a clear oscillatory 
behaviour of the populations of all four renormalized 
excitonic eigenstates, confirming that populations and 
coherences are coupled. 

Weak-Coupling Limit: In the limit that the system- 
bath coupling is sufficiently weak we may approximate 
the Markovian rates given in equation ( 45 ) by expand- 

of 



mg 



fC mn ,pq(t)- To first order, we find 



the correlation functions Cj£p (t) 



m powers 



-,(4) 



(W) 



,(3) 



where 



r iw) 



(mn) (pq) 



mnpq n 



^pq^mn^pq (^) 5 



(50) 



(51) 



is the weak-coupling bath correlation function. 

The rates in Eq. ( 45 ) can now be evaluated in 



straightforward manner by substituting in the form of the 



fcmn,pq(t) (as defined in Appendix |A[) into Eq. (51), and 
making use of the relation dte 1 = tt5(w) + iP[l/ui) 



fj, a ,an)Paa{t)??h.eT:e P denotes the principal value, to perform the in- 
tegrals over time. Hence, in the weak coupling limit, all 
four rates can be written in a generic form 



(48) 



Meanwhile, the off-diagonal matrix element p^ v (t) evolve 
as follows: 



ds 



r {w) , s 



"/a/3, [iv 

(e) - iS a p !fJ , v (e). 



(52) 
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(2) (3) 

only the rates T a ^ aj3 and T a p a p have non-zero contri- 
bution and evaluate to: 



r 



(2) 

a/3,a/3 



,(3) 



E V* n 5 ma 6 nfS T( s \e a0 ) 

(mn) 

E 

(mn) 



■ a p,a/3- / ^y^Jm a 5nflT {S \-e a p), (55) 

where we have defined the strong coupling limit rate as 
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FIG. 6. Dynamics of the eigenstate populations using the 
secular Markovian polaron master equation evaluated in the 
polaron frame (dash-dotted) and the lab frame (solid). 



It can be shown that in this limit, the site populations fol- 
low the Pauli master equation presented in Eq. ( 48 ) . On 



substituting in the forms of the rates T^ 2 j a/3 and T [ aj3 
from above, and after performing some simple manip- 
ulations, we arrive at the strong system-bath coupling 
master equation: 



"(3) 



dpnnjt) 

dt 



E 

m^n 

E 

m^n 



V r ^ m Re[T^(e nm )}p mm (t) 
V^ n Re[T^\e mn )]p nn (t). (57) 



Notice that the above is exactly the incoherent Forster 
dynamics with rates V^ n Re[T^ (e mn )] corresponding to 
the Forster transfer rate from site m to site n. 



V. CONCLUDING REMARKS 



Here 



laP,nv{ e ) — 2 E E MgBuv PmnPpq n ^™«,P<j( e ) 
(mn) (pq) 

x (coth (/3e/2) - 1) (53) 
is the expected single-phonon relaxation rate, while 

S a p,nv{t) = E E ^« B uv Pmnfipq 



{ran) {pq) 



XP 



J(w) A mn p q(w) , . . 

— ,, 2 _ g 2 (^-ecoth/3o;/2) 



(54) 



is the associated bath-induced energy shift. Hence, pro- 
vided that we may legitimately perform the expansion 
in /C mn>M (t), we see that our master equation should 
correctly capture the expected Redfield dynamics of the 
system in the weak-coupling limit. 

Forster Limit: In the opposite regime of strong system- 
bath interaction, the bath renormalisation factors j3 mn 
tend to zero. In this limit, electronic couplings are renor- 
malised to zero and the eigenstate basis is simply the site 
basis (u rna = S ma ). It can be shown that when fl mn — > 0, 



It is currently of much interest to develop theories 
of multichromopore electronic excitation dynamics ca- 
pable of bridging the gap between the limiting cases of 
weak and strong exciton-phonon coupling, whilst remain- 
ing computationally tractable as the number of chro- 
mophores increases. In this context, modified pertur- 
bative methodologies, as the one presented here, provide 
a valuable alternative to exact treatments. In particu- 
lar, this paper generalizes the polaron-modified pertur- 
bative master equation originally presented for a donor- 
aceptor pair to the case of multichromophore excita- 
tion dynamics. Explicit expressions for the homogeneous 
and inhomogeneous super-operators for an arbitrary non- 
equilibrium bath initial state have been presented. To 
illustrate the scope of this many-site theory, we have in- 
vestigated electronic excitation dynamics in a four-site 
subsystem of the FMO complex under the influence of a 
structured phonon bath that includes a localised vibra- 
tional mode. Our results indicate that in this example 
the non-equilibrium bath dynamics, captured by the in- 
homogeneous contribution, is crucial to give an accurate 
account of the origin and time scale of oscillations on the 
ultrafast scale. In particular, we show how the theory can 
describe the enhancement and modification of the oscilla- 
tory dynamics due to strong coupling to a localised, yet 
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broadened vibrational mode. In a separate publication 
we will be describing in detail how the interplay between 
electronic coherence and localised vibrational modes can 
generate rich behaviour of transfer of excitions in light 
harvesting systems. The ability to understand this inter- 
play will give insights into the role of coherent dynamics 
in systems where the excitonic transitions may be reso- 
nant with localised vibrational modes. 

In addition to calculating site population dynamics in 
the lab frame, we have outlined a framework for evaluat- 
ing all possible electronic observables in the lab frame. 
This will allow a full reconstruction of the lab frame 
density matrix for excitation dynamics and hence enable 
comparisons with experimental observations. However, 
explicit calculation of non-equilibrium contributions to 
the expected values of certain electronic operators are be- 
yond the scope of this paper, so we have presented here 
a zeroth-order approximation of such contributions and 
leave the full calculation for a forthcoming publication. 

Finally, we have presented the Markovian and secular 
approximations of the multichromophore polaron- 
modified master equation. Despite these approximations 
being satisfied in the polaron frame, we have shown that 
both approximations do not necessarily have to hold in 
the untransformed lab frame. This suggests that there 
is scope for Markov polaron master equations to capture 
aspects of non-Markovian or non-secular dynamics in 
the lab frame in particular parameter regimes. 

Note: on completation of this work we were made 
aware of a similar work by Jang 64 . Our approaches are 
similar in scope and objective but our analysis present 
complementary insights. 
Acknowledgements 

We would like to thank Michael Thorwart for discus- 
sions and comments to our paper and to Seogjoo Jang 
for providing a pre-print of his manuscript while ours 
was under revision. AK and AOC acknowledge funding 
from the EPSRC (grant No. EP/G005222/1 ). AN is 
grateful with the EPSRC and Imperial College London 
for support. 



Appendix A: Homogeneous correlations functions 

In order to calculate the bath correlation functions, 
the polaron frame bath operators are written in terms of 
displacement operators: 



B mn (t) = Y[D{5a ktmn (t)) 



(Al) 



where the displacement operator is in general defined as 
^(c^k) = e ak *~ a * k . Using the following properties of 
displacement operators 



D(a k )D(p k ) = e^-<^' 2 D{a k + /3 k ) 
(D(a k )) =exp(- i|a k | 2 coth(/3u; k /2)) (A2) 

the final expressions for the homogeneous bath correla- 
tion functions become 



(B mn (t)B pq ( S )) | = p mMe -K mn . pq{ t-s) lh 

{B} nn {t)B pq {s)) \ _ (t-s) _ ^ 



(Bl n (t)Bl q (s)) 
(B mn (t)Bl q (s)} 



(A3) 



where the correlation function JC mntPq (t) is defined as 



fc m n,pq(t) — j doj^-^-X mn<pq ( coth(#<;/2) cos(w(t)) 



—i sin(u>(t)) , 

(A4) 

and the spatial correlation function is defined as \ m n, P q = 
A TOiP - A m>g ~ A n<p + A„ !? . Here A TOiP describes the 
degree of spatial correlation between sites m and p. 
For the propagating modes model of spatial correlations 
A TOiP = A mjP (cj) = sm.c(u)d mn /v p h). 



Appendix B: Inhomogeneous correlation functions 

Let us consider a general initial state within the lab 
frame (i.e. prior to polaron tranformation): x(0) = 
Eij Pij{Q) a i a j ®Pb, where p B denotes the thermal equi- 
librium bath state in the lab frame. Transforming into 
the polaron frame we find the initial state 

X(0) =Y / p ij (0)a+aT]l^D(a k , i )p B D(-a k , j ). 

ij k 

(Bl) 

Here py(0) = (3ijpij(Q) and denotes the ij'th clement of 
the initial system density operator in the polaron frame. 

The irrelevant part of the total system-bath density 
matrix at time zero is then given by 



(°)^ +(T 7 II (^D{a k ^p B D{-a ktj ) - p B 



ij k 

Pij (°) °t <*j QijPB- (BS 



Notice that we have defined QijpB as the state account- 
ing for the difference between the displaced bath and the 
bath thermal equilibrium in the polaron frame. 

Expectation values with respect to the state QijPs can 
be expressed in terms of expectation values taken with 



14 



respect to the thermal equilibrium state in the polaron 
frame, as follows: 

( x )q»pb = l[(D(-a Kj )XD(a kti ))- PB - (X) PB (B3) 

k 

Using this identity and the previous properties of dis- 
placement operators, we can now calculate the various 
inhomogeneous correlation functions. The full expres- 
sion for the correlation function appearing in the first 
order term of the inhomogeneous super-operator can be 
evaluated as: 

(BmnWJQijPB = Pmn(fij,mn\t) ~ !)• (B4) 

Here the correlation function fij im .n{t) is defined as 



f (+\ — ~ Io° ^T^^j.in" coth(^w/2) cos(uit) 



xe 



sin(wi) 



(B5) 



The spatial correlation factor A.y ;mn is as defined in Ap- 
pendix [XJ while a second spatial correlation function is 



introduced: A: 

The correlation functions appearing in the second order 
inhomogeneous term are given by: 



{B mn {t)B pq (s)) Q i]jiB 

= 8 8 ((f (t)f- (s) - iV^^piC- 8 ) 

fij,mn(t) fij,pqi s ) ^ 5 

{Bl n (t)B pq {s)) Qi] p B 

= f3rnn(3pq(^(fij, mn (t)fij, P q(s) — l)e mn ' pq ( ) 

{Bmn{t)Bl q {s))Q i]i>B 

— finmfipq (^(fij,mn(t) fij^ pq (s) — l)e mn - p9 ^ ) 

~fij,mn{t) ~~ fij,pq( s ) + 2^ , 

(Bl n (t)Bl q (s)) QijPB 

— Prnn8 p q(^{fij^ mn (t)f i j pq (s)~l^e mn.pqi ) 
~ fij,mn(t) ~~ fij,pqi s ) • 



(B6) 



Here we have introduced a final correlation function: 

fij,mn(t) (fij,mn{t)) 

_ /" du^^Xi^rnn coth(/3w/2) cos(wt) 



Appendix C: Numerical Integration 

For convenience we numerically solve the dynamics 
in the polaron frame within the eigenstate basis of the 
renormalised Hamiltonian in Equation [4j In this basis, 
we may write the polaron master equation as: 



- — ^2 W/V (*) + lap (*) 



(CI) 



Here, R a p,ii.v{t) and I a p{t) are time-dependent tensors 
corresponding to the homogeneous and inhomogeneous 
superoperators respectvely. To simplify the numerics, we 
flatten the system density matrix to form a vector de- 
scribing the state: p = {fin, P12, P13, ■ ■ ■ , Pnn) T ■ In this 
new representation, we may write the master equation in 
terms of the following matrix equation: 



-p(t) = R{t).p(t) 



Lit) 



(C2) 



x e 



This matrix equation is numerically integrated using the 
fourth order Runge-Kutta method. 

At each time step during the numerical integration, 
the elements of the homogeneous matrix R(t) and the 

inhomogeneous vector I(t) are determined from the ex- 
pressions in Equations [l8j [24] and 26 The most compu- 
tationally intensive step in evaluating these two terms 
occurs in performing the integrations to calculate the 
time-dependent homogeneous and inhomogeneous rates. 
To reduce operation time, at each time step we calculate 
all rates first, before building up the homogeneous ma- 
trix and inhomogeneous vector. Furthermore, we notice 
that at each time step all the rates can all be calculated 
independently. Therefore, we may also utilise paralleli- 
sation algorithms to further enhance the performance of 
the numerical integration. 
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